0. Setup

Load the data from Study 2.

Each row is a different participant in Study 2 (there were 80). Columns 1-68 are variables that summarise the whole session. Columns 69-493 refer to events in each round of the Visible condition. Columns 494-918 refer to the Hidden condition. This experiment followed a within-subjects design (i.e. participants completed both the Hidden and Visible conditions). In the counterbalancing column, 0 means that the participant played the Visible game before the Hidden game, and 1 means that the participant played the Hidden game before the Visible game.

1. Binomial - Number of shocks

Trim the dataset and get it in long-format.

As in Study 1, we first test whether the probability of a shock occurring differs between the two conditions. Although this experiment followed a within-subjects design, random effects are not necessary here because the occurrence of shocks are independent between the two conditions.

All analyses in this document will use the brm() function. Learn more about the brms package here.

## function(d2.01) {
##   brm(data = d2.01, family = binomial,
##       overall_num_shocks | trials(rounds_survived) ~ 0 + Intercept + Condition,
##       prior = c(prior(normal(0, 1), class = b)),
##       iter = 2000, warmup = 1000, chains = 4, cores = 4,
##       seed = 2113)
## }

We set the seed to a random number, to make the results reproducible. Here are the priors that were set for this model.

Let’s look at the results.

##  Family: binomial 
##   Links: mu = logit 
## Formula: overall_num_shocks | trials(rounds_survived) ~ 0 + Intercept + Condition 
##    Data: d2.01 (Number of observations: 160) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -1.36      0.06    -1.47    -1.24 1.00     1562     1870
## Condition    -0.04      0.09    -0.21     0.13 1.00     1511     1854
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plotting the parameters.

The Condition parameter is -0.04, and its 95% credible intervals cross 0, implying that the probability of a shock is no different across the two conditions. On the probability scale:

post <- readd(post2.01)

visible_prob <- inv_logit_scaled(post$b_Intercept)
hidden_prob  <- inv_logit_scaled(post$b_Intercept + post$b_Condition)
difference   <- hidden_prob - visible_prob

quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
## -0.03 -0.01  0.02

2. Gaussian - Total amount of cattle lost due to shocks

Trim the dataset again.

We now fit a model to determine if the total amount of cattle lost due to shocks varies between conditions. Again, no random effects are needed because the outcome is independent across conditions (generated stochastically by the game), despite the within-subject design.

## function(d2.02) {
##   brm(data = d2.02, family = gaussian,
##       total_cattle_lost ~ 0 + Intercept + Condition,
##       prior = c(prior(normal(0, 100), class = b, coef = 'Intercept'),
##                 prior(normal(0, 5), class = b, coef = 'Condition')),
##       iter = 2000, warmup = 1000, chains = 4, cores = 4,
##       seed = 2113)
## }

The priors we used.

The results.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: total_cattle_lost ~ 0 + Intercept + Condition 
##    Data: d2.02 (Number of observations: 160) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    51.29      2.13    47.09    55.47 1.00     2258     2302
## Condition    -6.79      2.77   -12.16    -1.35 1.00     2269     2347
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    20.42      1.17    18.31    22.84 1.00     2446     2197
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plot the parameters.

On average, 7 more cattle were lost in the Hidden condition than the Visible condition.

3. Binomial - Probability of requesting

Get a long-format data frame with binary request decisions over all rounds, for both conditions. If request == NA, player has died and been removed from the game, so we drop those rows.

This leaves us with 2834 request decisions.

We now fit a varying intercept and slope model, with participants nested within groups. We allow the slopes for both round number and condition to vary, to fit with the experiment’s within-subjects design (participants completed multiple rounds, in both conditions).

## function(d2.03) {
##   brm(data = d2.03, family = bernoulli,
##       request ~ 0 + Intercept + round_number + Condition
##       + (0 + Intercept + round_number + Condition | Group/ID),
##       prior = c(prior(normal(0, 1), class = b),
##                 prior(student_t(3, 0, 10), class = sd),
##                 prior(lkj(1), class = cor)),
##       iter = 2000, warmup = 1000, chains = 4, cores = 4,
##       sample_prior = TRUE,
##       control = list(adapt_delta = 0.99),
##       seed = 2113)
## }

Here are the priors for the model we just fitted.

Now let’s see the results.

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: request ~ 0 + Intercept + round_number + Condition + (0 + Intercept + round_number + Condition | Group/ID) 
##    Data: d2.03 (Number of observations: 2834) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   0.63      0.28     0.07     1.15 1.01      448      807
## sd(round_number)                0.05      0.02     0.01     0.08 1.01      388      408
## sd(Condition)                   1.62      0.32     1.01     2.29 1.00      993     1153
## cor(Intercept,round_number)    -0.37      0.40    -0.88     0.63 1.01      499      940
## cor(Intercept,Condition)       -0.24      0.34    -0.80     0.56 1.01      373      361
## cor(round_number,Condition)     0.35      0.30    -0.31     0.87 1.02      312      396
## 
## ~Group:ID (Number of levels: 80) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   0.71      0.25     0.22     1.18 1.01      376      714
## sd(round_number)                0.04      0.02     0.00     0.08 1.01      271      980
## sd(Condition)                   1.05      0.31     0.50     1.71 1.01      644     1328
## cor(Intercept,round_number)    -0.44      0.42    -0.91     0.66 1.00      624     1242
## cor(Intercept,Condition)       -0.33      0.32    -0.82     0.41 1.01      500     1238
## cor(round_number,Condition)     0.16      0.40    -0.67     0.85 1.01      422      960
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       -1.00      0.18    -1.37    -0.66 1.00     3278     2683
## round_number    -0.04      0.01    -0.07    -0.02 1.00     2897     2569
## Condition        0.63      0.31     0.03     1.22 1.00     2935     2619
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plotting the parameters.

Finally, let’s see the trace plots.

Looks like Stan sampled efficiently.

In this model, the fixed effect of condition is 0.63, with 95% credible intervals above 0. This implies that the participants are more likely to request from their partner in the Hidden condition. Converting to the probability scale:

post <- readd(post2.03)

visible_prob <- inv_logit_scaled(post$b_Intercept)

visible_prob %>%
  median() %>%
  round(2)
## [1] 0.27
hidden_prob  <- inv_logit_scaled(post$b_Intercept + post$b_Condition)

hidden_prob %>%
  median() %>%
  round(2)
## [1] 0.41
difference <- hidden_prob - visible_prob

quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
##  0.01  0.14  0.28

The absolute probability difference between the conditions is +0.14 (median), with 95% CIs above 0. Participants were more likely to request from their partner in the Hidden condition.

We compute a Bayes factor for this difference between probabilities. This is the inverse of the alternative hypothesis that the two conditions are equal.

## [1] 2.84

This Bayes factor implies weak support for the hypothesis that the probabilities differ across conditions.

4a. Binomial - Probability of requesting when above the minimum threshold

Get a long-format data frame with binary request decisions over all rounds, for both conditions. If request == NA, player has died and been removed from the game, so we drop those rows. However, we also filter out rows in which the player was below the minimum survival threshold (64 cattle).

This leaves us with 2368 request decisions.

We now fit a varying intercept and slope model, with participants nested within groups. Again, we allow the slopes for both round number and condition to vary.

## function(d2.04) {
##   brm(data = d2.04, family = bernoulli,
##     request ~ 0 + Intercept + round_number + Condition
##     + (0 + Intercept + round_number + Condition | Group/ID),
##     prior = c(prior(normal(0, 1), class = b),
##               prior(student_t(3, 0, 10), class = sd),
##               prior(lkj(1), class = cor)),
##     iter = 2000, warmup = 1000, chains = 4, cores = 4,
##     sample_prior = TRUE,
##     control = list(adapt_delta = 0.99),
##     seed = 2113)
## }

The priors.

The results.

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: request ~ 0 + Intercept + round_number + Condition + (0 + Intercept + round_number + Condition | Group/ID) 
##    Data: d2.04 (Number of observations: 2368) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   1.12      0.47     0.16     2.03 1.00      548      707
## sd(round_number)                0.08      0.03     0.01     0.14 1.01      460      777
## sd(Condition)                   1.93      0.75     0.27     3.33 1.01      548      889
## cor(Intercept,round_number)    -0.14      0.43    -0.81     0.77 1.01      600     1478
## cor(Intercept,Condition)       -0.04      0.41    -0.75     0.77 1.00      674     1477
## cor(round_number,Condition)     0.33      0.38    -0.60     0.89 1.00      648      870
## 
## ~Group:ID (Number of levels: 80) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   1.54      0.36     0.88     2.28 1.00      718     1407
## sd(round_number)                0.09      0.03     0.04     0.14 1.00      546      901
## sd(Condition)                   2.29      0.64     1.19     3.65 1.01      663     1381
## cor(Intercept,round_number)    -0.60      0.24    -0.92     0.03 1.01      947     1330
## cor(Intercept,Condition)       -0.26      0.26    -0.71     0.27 1.00      935     1533
## cor(round_number,Condition)     0.51      0.29    -0.17     0.93 1.00      444      982
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       -1.98      0.32    -2.64    -1.38 1.00     2653     2709
## round_number    -0.10      0.03    -0.16    -0.05 1.00     2649     2713
## Condition        0.19      0.49    -0.83     1.09 1.00     3187     3003
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Creating a forest plot of parameters.

Trace plots.

Looks like Stan sampled efficiently.

In this model, the fixed effect of condition is 0.19, with 95% credible intervals crossing 0. This implies that, at least when above the minimum survival threshold, the average probability of requesting did not differ between conditions.

Converting the fixed effects to the probability scale:

post <- readd(post2.04a)

visible_prob <- inv_logit_scaled(post$b_Intercept)

visible_prob %>%
  median() %>%
  round(2)
## [1] 0.12
hidden_prob  <- inv_logit_scaled(post$b_Intercept + post$b_Condition)

hidden_prob %>%
  median() %>%
  round(2)
## [1] 0.14
difference <- hidden_prob - visible_prob

quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
## -0.08  0.02  0.17

The absolute probability difference between the conditions is +0.02 (median), with 95% CIs above 0. Participants were no more likely to request from their partner in the Hidden condition.

We compute a Bayes factor for this difference between probabilities.

## [1] 0.32

This Bayes factor implies moderate support for the hypothesis that the probabilities are equal in each condition.

4b. Binomial - Probability of requesting when above the minimum threshold (interacting with counterbalancing)

We found no effect of condition in the previous model. Is this because of order effects?

## function(d2.04) {
##   brm(data = d2.04, family = bernoulli,
##       request ~ 0 + Intercept + round_number + Condition*Counterbalancing
##       + (0 + Intercept + round_number + Condition | Group/ID),
##       prior = c(prior(normal(0, 1), class = b),
##                 prior(student_t(3, 0, 10), class = sd),
##                 prior(lkj(1), class = cor)),
##       iter = 2000, warmup = 1000, chains = 4, cores = 4,
##       sample_prior = TRUE,
##       control = list(adapt_delta = 0.99),
##       seed = 2113)
## }

The priors.

The results.

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: request ~ 0 + Intercept + round_number + Condition * Counterbalancing + (0 + Intercept + round_number + Condition | Group/ID) 
##    Data: d2.04 (Number of observations: 2368) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   1.13      0.49     0.13     2.06 1.01      620      682
## sd(round_number)                0.08      0.03     0.01     0.14 1.00      618      779
## sd(Condition)                   1.30      0.70     0.09     2.72 1.00      591     1160
## cor(Intercept,round_number)    -0.13      0.42    -0.80     0.80 1.00      863     1499
## cor(Intercept,Condition)       -0.15      0.45    -0.86     0.80 1.00      878     1524
## cor(round_number,Condition)     0.26      0.43    -0.72     0.90 1.01      861     1926
## 
## ~Group:ID (Number of levels: 80) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   1.53      0.35     0.86     2.25 1.00      964     2090
## sd(round_number)                0.08      0.02     0.04     0.14 1.00      882     1531
## sd(Condition)                   2.16      0.53     1.20     3.28 1.00      846     1902
## cor(Intercept,round_number)    -0.62      0.23    -0.93    -0.05 1.00     1411     2184
## cor(Intercept,Condition)       -0.30      0.24    -0.71     0.21 1.00     1351     2267
## cor(round_number,Condition)     0.55      0.26    -0.05     0.94 1.00      571     1184
## 
## Population-Level Effects: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept                     -2.09      0.40    -2.88    -1.32 1.00     3386     3295
## round_number                  -0.10      0.03    -0.15    -0.05 1.00     2992     2940
## Condition                     -0.38      0.50    -1.38     0.58 1.00     4079     3339
## Counterbalancing               0.21      0.48    -0.76     1.12 1.00     3580     3035
## Condition:Counterbalancing     1.55      0.65     0.24     2.76 1.00     3334     2696
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Creating a forest plot of parameters.

Trace plots.

Looks like Stan sampled efficiently.

In this model, the interaction parameter is 1.55, with 95% CIs above zero. This implies that there is an interaction effect. Condition only has an effect on cheating when participants played the Hidden game first.

Converting the fixed effects to the probability scale:

post <- readd(post2.04b)

visible_VisibleFirst_prob <- inv_logit_scaled(post$b_Intercept)

visible_VisibleFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.11
hidden_VisibleFirst_prob  <- inv_logit_scaled(post$b_Intercept + post$b_Condition)

hidden_VisibleFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.08
difference_VisibleFirst <- hidden_VisibleFirst_prob - visible_VisibleFirst_prob

quantile(difference_VisibleFirst,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
## -0.11 -0.03  0.07

When participants play the Visible game first, the absolute probability difference between the conditions is -0.03 (median), with 95% CIs crossing 0. Participants were no more likely to request from their partner in the Hidden condition.

We compute a Bayes factor for this difference between probabilities.

## [1] 0.31

This Bayes factor implies moderate support for the hypothesis that the probabilities are equal in each condition.

What about for when participants play the Hidden game first?

visible_HiddenFirst_prob <- inv_logit_scaled(post$b_Intercept + post$b_Counterbalancing)

visible_HiddenFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.13
hidden_HiddenFirst_prob  <- inv_logit_scaled(post$b_Intercept + post$b_Counterbalancing + 
                                          post$b_Condition + post$`b_Condition:Counterbalancing`)

hidden_HiddenFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.34
difference_HiddenFirst <- hidden_HiddenFirst_prob - visible_HiddenFirst_prob

quantile(difference_HiddenFirst,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
##  0.00  0.20  0.43

When participants play the Hidden game first, the absolute probability difference between the conditions is +0.20 (median), with 95% CIs above 0. Participants are more likely to request from their partner when beneath the threshold in the Hidden condition.

We compute a Bayes factor for this difference between probabilities.

## [1] 2.44

This Bayes factor implies weak support for the hypothesis that the probabilities differ across conditions.

5. Binomial - Probability of not responding to a request

Get long-format data frame with ‘received’ variable (i.e. how much a player received on any given round). We swap this around so it reflects how much the player gave to their partner (i.e. how much their partner received). We drop rows with NAs, since partners did not request help in that particular round. We then code whether the player gave nothing in response to the request (1) or gave at least one cattle (0).

This leaves us with 716 possible responses to requests.

We then fit the varying intercept and slope model, again with participants nested within groups.

## function(d2.05) {
##   brm(data = d2.05, family = bernoulli,
##       notResponded ~ 0 + Intercept + round_number + Condition
##       + (0 + Intercept + round_number + Condition | Group/ID),
##       prior = c(prior(normal(0, 1), class = b),
##                 prior(student_t(3, 0, 10), class = sd),
##                 prior(lkj(1), class = cor)),
##       iter = 2500, warmup = 1000, chains = 4, cores = 4,
##       sample_prior = TRUE,
##       control = list(adapt_delta = 0.99),
##       seed = 2113)
## }

The priors we used.

The results.

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: notResponded ~ 0 + Intercept + round_number + Condition + (0 + Intercept + round_number + Condition | Group/ID) 
##    Data: d2.05 (Number of observations: 716) 
## Samples: 4 chains, each with iter = 2500; warmup = 1000; thin = 1;
##          total post-warmup samples = 6000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   1.04      0.48     0.11     2.03 1.00     1087      897
## sd(round_number)                0.09      0.04     0.01     0.18 1.00      649      871
## sd(Condition)                   1.19      0.56     0.16     2.33 1.00      538      853
## cor(Intercept,round_number)     0.09      0.42    -0.68     0.85 1.00      683     1140
## cor(Intercept,Condition)       -0.40      0.41    -0.92     0.62 1.00      709     1490
## cor(round_number,Condition)     0.37      0.41    -0.60     0.93 1.00      991     1565
## 
## ~Group:ID (Number of levels: 78) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   0.52      0.41     0.02     1.49 1.00      664     1547
## sd(round_number)                0.03      0.02     0.00     0.09 1.00     1961     2966
## sd(Condition)                   0.55      0.45     0.02     1.65 1.00      595     1337
## cor(Intercept,round_number)    -0.09      0.51    -0.90     0.87 1.00     2710     3871
## cor(Intercept,Condition)       -0.41      0.50    -0.98     0.74 1.00      987     2731
## cor(round_number,Condition)    -0.05      0.50    -0.90     0.87 1.00     4005     4462
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       -1.70      0.37    -2.49    -1.01 1.00     2261     3752
## round_number    -0.06      0.04    -0.14    -0.00 1.00     1347     2773
## Condition        0.49      0.38    -0.28     1.21 1.00     2662     3868
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plotting the parameters.

Let’s look at the trace plots to make sure Stan sampled efficiently.

Rhat values, n_eff values, and trace plots look okay.

The fixed effect of condition is 0.49, with 95% credible intervals that cross 0. This implies that participants were no more likely to not respond to requests in the Hidden condition.

As before, we sample from the posterior and convert this difference to the absolute probability scale.

post <- readd(post2.05)

visible_prob <- inv_logit_scaled(post$b_Intercept)

visible_prob %>%
  median() %>%
  round(2)
## [1] 0.16
hidden_prob  <- inv_logit_scaled(post$b_Intercept + post$b_Condition)

hidden_prob %>%
  median() %>%
  round(2)
## [1] 0.23
difference <- hidden_prob - visible_prob

quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
## -0.04  0.07  0.19

The probability difference between the conditions is 0.07 (median), with 95% CIs crossing 0.

We compute a Bayes factor for this difference between probabilities.

## [1] 0.75

This Bayes factor implies anecdotal support for the hypothesis that the probabilities are equal across conditions.

6a. Binomial - Probability of fulfilling a request when able

Again, the data wrangling for this model is a little trickier.

  1. First, we get a long-format data frame with (a) the player’s herd size that round, (b) how much the player received that round, and (c) how much the player requested that round.
  2. Next, we flip the latter two variables to get (b’) how much the player gave, and (c’) how much the player’s partner requested. ‘Flipping’ is possible because partners sit next to the focal player in the data frame.
  3. We drop rows with NAs, as no requesting happened this round.
  4. We keep only rows in which the partner’s request could be fulfilled without dropping the player below the minimum survival threshold (i.e. the player was able to give).
  5. Finally, we code whether the player fulfilled the request by giving what was asked or more (0) or did not fulfill the request (1).

This leaves us with 473 possible response decisions in which the player was able to give their partner what they asked for. Our outcome variable is whether they fulfilled that request or not.

We now fit the varying intercept and slope model, again with participants nested within groups.

## function(d2.06) {
##   brm(data = d2.06, family = bernoulli,
##       notFulfilled ~ 0 + Intercept + round_number + Condition
##       + (0 + Intercept + round_number + Condition | Group/ID),
##       prior = c(prior(normal(0, 1), class = b),
##                 prior(student_t(3, 0, 10), class = sd),
##                 prior(lkj(1), class = cor)),
##       control = list(adapt_delta = 0.99),
##       sample_prior = TRUE,
##       iter = 2500, warmup = 1000, chains = 4, cores = 4,
##       seed = 2113)
## }

Here are the priors we used.

Let’s see the results.

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: notFulfilled ~ 0 + Intercept + round_number + Condition + (0 + Intercept + round_number + Condition | Group/ID) 
##    Data: d2.06 (Number of observations: 473) 
## Samples: 4 chains, each with iter = 2500; warmup = 1000; thin = 1;
##          total post-warmup samples = 6000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   1.22      0.56     0.15     2.38 1.00     1367     1065
## sd(round_number)                0.16      0.07     0.04     0.33 1.00      915     1261
## sd(Condition)                   0.81      0.61     0.03     2.27 1.00     1813     2487
## cor(Intercept,round_number)     0.21      0.40    -0.57     0.90 1.00     1394     1378
## cor(Intercept,Condition)       -0.05      0.49    -0.88     0.86 1.00     3517     3473
## cor(round_number,Condition)     0.22      0.47    -0.76     0.93 1.00     3569     4125
## 
## ~Group:ID (Number of levels: 73) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   0.64      0.47     0.02     1.77 1.00     1628     2715
## sd(round_number)                0.08      0.05     0.01     0.19 1.00     1509     2655
## sd(Condition)                   0.77      0.56     0.03     2.09 1.00     1428     2773
## cor(Intercept,round_number)    -0.19      0.49    -0.92     0.82 1.00     2488     3703
## cor(Intercept,Condition)       -0.30      0.49    -0.95     0.77 1.00     2710     3777
## cor(round_number,Condition)     0.07      0.50    -0.86     0.89 1.00     2855     4481
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       -0.54      0.39    -1.34     0.19 1.00     3069     4056
## round_number    -0.13      0.06    -0.26    -0.04 1.00     2610     3225
## Condition        0.94      0.40     0.12     1.68 1.00     4551     4170
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plotting the parameters.

Let’s look at the trace plots to make sure Stan sampled efficiently.

Stan sampled efficiently.

The fixed effect of condition is 0.94, with 95% credible intervals above 0. This implies that participants were more likely to not fulfill requests (when able) in the Hidden condition.

On the absolute probability scale.

post <- readd(post2.06a)

visible_prob <- inv_logit_scaled(post$b_Intercept)

visible_prob %>%
  median() %>%
  round(2)
## [1] 0.37
hidden_prob  <- inv_logit_scaled(post$b_Intercept + post$b_Condition)

hidden_prob %>%
  median() %>%
  round(2)
## [1] 0.6
difference <- hidden_prob - visible_prob

quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
##  0.03  0.23  0.39

The probability difference between the conditions is 0.23 (median), with 95% CIs above zero.

We compute a Bayes factor for this difference between probabilities.

## [1] 6.44

This Bayes factor implies moderate support for the hypothesis that the probabilities differ across conditions.

6b. Binomial - Probability of fulfilling request when able (interacting with counterbalancing)

As before, we test for order effects by including an interaction.

## function(d2.06) {
##   brm(data = d2.06, family = bernoulli,
##       notFulfilled ~ 0 + Intercept + round_number + Condition*Counterbalancing
##       + (0 + Intercept + round_number + Condition | Group/ID),
##       prior = c(prior(normal(0, 1), class = b),
##                 prior(student_t(3, 0, 10), class = sd),
##                 prior(lkj(1), class = cor)),
##       control = list(adapt_delta = 0.99),
##       sample_prior = TRUE,
##       iter = 2500, warmup = 1000, chains = 4, cores = 4,
##       seed = 2113)
## }

Here are the priors we used.

Let’s see the results.

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: notFulfilled ~ 0 + Intercept + round_number + Condition * Counterbalancing + (0 + Intercept + round_number + Condition | Group/ID) 
##    Data: d2.06 (Number of observations: 473) 
## Samples: 4 chains, each with iter = 2500; warmup = 1000; thin = 1;
##          total post-warmup samples = 6000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   1.17      0.56     0.15     2.37 1.01     1375     1383
## sd(round_number)                0.17      0.07     0.04     0.33 1.01     1037     1258
## sd(Condition)                   0.83      0.63     0.04     2.33 1.00     1061     2277
## cor(Intercept,round_number)     0.22      0.38    -0.53     0.88 1.00     1118     2273
## cor(Intercept,Condition)       -0.09      0.49    -0.88     0.86 1.00     2739     3326
## cor(round_number,Condition)     0.23      0.46    -0.77     0.91 1.00     2979     4584
## 
## ~Group:ID (Number of levels: 73) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                   0.70      0.49     0.03     1.86 1.00     1062     1754
## sd(round_number)                0.08      0.05     0.00     0.20 1.00     1089     2000
## sd(Condition)                   0.86      0.61     0.04     2.29 1.00     1190     2541
## cor(Intercept,round_number)    -0.19      0.50    -0.92     0.82 1.00     1819     3484
## cor(Intercept,Condition)       -0.32      0.49    -0.95     0.77 1.00     1889     3232
## cor(round_number,Condition)     0.12      0.50    -0.84     0.90 1.00     2692     3945
## 
## Population-Level Effects: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept                     -0.76      0.44    -1.67     0.08 1.00     4018     4474
## round_number                  -0.13      0.05    -0.25    -0.04 1.00     2512     3362
## Condition                      0.63      0.51    -0.42     1.58 1.00     4191     4119
## Counterbalancing               0.44      0.55    -0.67     1.49 1.00     4565     4552
## Condition:Counterbalancing     0.48      0.59    -0.67     1.66 1.00     5237     4527
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plotting the parameters.

The interaction parameter has 95% CIs that cross zero, indicating no interaction effect.

Converting the fixed effects to the probability scale:

post <- readd(post2.06b)

visible_VisibleFirst_prob <- inv_logit_scaled(post$b_Intercept)

visible_VisibleFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.32
hidden_VisibleFirst_prob  <- inv_logit_scaled(post$b_Intercept + post$b_Condition)

hidden_VisibleFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.47
difference_VisibleFirst <- hidden_VisibleFirst_prob - visible_VisibleFirst_prob

quantile(difference_VisibleFirst,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
## -0.09  0.15  0.36

When participants play the Visible game first, the absolute probability difference between the conditions is +0.15 (median), with 95% CIs crossing 0. Participants were no more likely to not fulfill their partner’s request in the Hidden condition.

We compute a Bayes factor for this difference between probabilities.

## [1] 1.38

This Bayes factor implies anecdotal support for the hypothesis that the probabilities are different across conditions.

What about for when participants play the Hidden game first?

visible_HiddenFirst_prob <- inv_logit_scaled(post$b_Intercept + post$b_Counterbalancing)

visible_HiddenFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.43
hidden_HiddenFirst_prob  <- inv_logit_scaled(post$b_Intercept + post$b_Counterbalancing + 
                                          post$b_Condition + post$`b_Condition:Counterbalancing`)

hidden_HiddenFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.69
difference_HiddenFirst <- hidden_HiddenFirst_prob - visible_HiddenFirst_prob

quantile(difference_HiddenFirst,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
##  0.03  0.26  0.46

When participants play the Hidden game first, the absolute probability difference between the conditions is +0.26 (median), with 95% CIs above 0. Participants are more likely to not fulfill their partner’s request in the Hidden condition.

We compute a Bayes factor for this difference between probabilities.

## [1] 6.88

This Bayes factor implies moderate support for the hypothesis that the probabilities differ across conditions.

7. Survival analysis - Condition

Get the data.

Fit the model.

## function(d2.07) {
##   brm(rounds_survived | cens(censored) ~ 0 + Intercept + Condition + (1 | Group/ID), 
##       data = d2.07, family = weibull, inits = "0",
##       prior = prior(normal(0, 2), class = b, coef = "Condition"),
##       iter = 2000, warmup = 1000, chains = 4, cores = 4,
##       control = list(adapt_delta = 0.9, max_treedepth = 15),
##       sample_prior = TRUE, seed = 2113)
## }

What priors did we use?

Results.

##  Family: weibull 
##   Links: mu = log; shape = identity 
## Formula: rounds_survived | cens(censored) ~ 0 + Intercept + Condition + (1 | Group/ID) 
##    Data: d2.07 (Number of observations: 160) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.39      0.24     0.02     0.92 1.00     1027     1515
## 
## ~Group:ID (Number of levels: 80) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.28      0.21     0.01     0.76 1.00     1178     1768
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     4.62      0.35     4.04     5.40 1.00     3271     2223
## Condition    -0.82      0.30    -1.43    -0.26 1.00     4174     2496
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     1.06      0.14     0.80     1.35 1.00     3160     2691
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plot parameters.

Trace plots.

Conditional effects.

Estimate survival rate in Visible condition.

post <- readd(post2.07)

set.seed(2113)
srateVisible <- rweibull(4000, post$shape, exp(post$b_Intercept))
median(srateVisible)
## [1] 65.89269

Estimate survival rate in Hidden condition.

set.seed(2113)
srateHidden <- rweibull(4000, post$shape, exp(post$b_Intercept + post$b_Condition))
median(srateHidden)
## [1] 29.91275
difference <- srateHidden - srateVisible
quantile(difference, c(0.025, 0.5, 0.975))
##        2.5%         50%       97.5% 
## -293.795972  -33.721701   -1.020057

Does survival rate differ across conditions? Get Bayes factor.

## [1] 5.93

Yes, it does.

8. Survival analysis - Proportion of Rule 1 cheating

Load data.

Fit the model.

## function(d2.08) {
##   brm(rounds_survived | cens(censored) ~ 0 + Intercept + prop_rule1 + (1 | Group/ID), 
##       data = d2.08, family = weibull, inits = "0",
##       prior = prior(normal(0, 2), class = b, coef = "prop_rule1"),
##       iter = 2000, warmup = 1000, chains = 4, cores = 4,
##       control = list(adapt_delta = 0.99),
##       sample_prior = TRUE, seed = 2113)
## }

What priors did we use?

Results.

##  Family: weibull 
##   Links: mu = log; shape = identity 
## Formula: rounds_survived | cens(censored) ~ 0 + Intercept + prop_rule1 + (1 | Group/ID) 
##    Data: d2.08 (Number of observations: 154) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.29      0.20     0.01     0.73 1.00     1700     2059
## 
## ~Group:ID (Number of levels: 80) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.23      0.17     0.01     0.62 1.00     1997     2183
## 
## Population-Level Effects: 
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      4.10      0.27     3.66     4.69 1.00     3870     2832
## prop_rule1     0.49      0.55    -0.52     1.66 1.00     6158     2010
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     1.13      0.16     0.84     1.45 1.00     3895     2973
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plot parameters.

Trace plots.

Conditional effects.

Does survival rate differ the more people cheat rule 1? Get Bayes factor.

## [1] 0.4

No. There is more evidence for the null hypothesis that the survival rate does not differ the more people break rule 1.

9. Survival analysis - Proportion of Rule 2 cheating

Load data.

Fit the model.

## function(d2.09) {
##   brm(rounds_survived | cens(censored) ~ 0 + Intercept + prop_rule2 + (1 | Group/ID), 
##       data = d2.09, family = weibull, inits = "0",
##       prior = prior(normal(0, 2), class = b, coef = "prop_rule2"),
##       iter = 2000, warmup = 1000, chains = 4, cores = 4,
##       control = list(adapt_delta = 0.99),
##       sample_prior = TRUE, seed = 2113)
## }

What priors did we use?

Results.

##  Family: weibull 
##   Links: mu = log; shape = identity 
## Formula: rounds_survived | cens(censored) ~ 0 + Intercept + prop_rule2 + (1 | Group/ID) 
##    Data: d2.09 (Number of observations: 117) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.28      0.21     0.01     0.80 1.00     1722     2106
## 
## ~Group:ID (Number of levels: 73) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.28      0.21     0.01     0.78 1.00     1395     1864
## 
## Population-Level Effects: 
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      4.26      0.36     3.71     5.09 1.00     3857     2367
## prop_rule2    -0.18      0.41    -0.99     0.64 1.00     4885     2980
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     1.39      0.26     0.93     1.95 1.00     4738     2646
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plot parameters.

Trace plots.

Conditional effects.

Does survival rate differ the more people cheat rule 2? Get Bayes factor.

## [1] 0.21

With a Bayes factor below 0.33, there’s strong support for the survival rate remaining the same as breaking rule 2 increases.

Session Info

sessionInfo()
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows >= 8 x64 (build 9200)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_New Zealand.1252  LC_CTYPE=English_New Zealand.1252    LC_MONETARY=English_New Zealand.1252
## [4] LC_NUMERIC=C                         LC_TIME=English_New Zealand.1252    
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tidyselect_1.0.0  forcats_0.4.0     stringr_1.4.0     dplyr_0.8.4       purrr_0.3.3       readr_1.3.1      
##  [7] tidyr_1.0.2       tibble_2.1.3      tidyverse_1.3.0   png_0.1-7         ggplot2_3.2.1     drake_7.10.0.9000
## [13] cowplot_1.0.0     brms_2.11.6       Rcpp_1.0.3        bayesplot_1.7.1  
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_1.4-1     ggridges_0.5.2       rsconnect_0.8.16     markdown_1.1         base64enc_0.1-3     
##   [6] fs_1.3.1             rstudioapi_0.11      farver_2.0.3         rstan_2.19.2         DT_0.12             
##  [11] fansi_0.4.1          mvtnorm_1.0-11       lubridate_1.7.4      xml2_1.2.2           bridgesampling_0.8-1
##  [16] knitr_1.28           shinythemes_1.1.2    jsonlite_1.6.1       broom_0.5.4          dbplyr_1.4.2        
##  [21] shiny_1.4.0          compiler_3.6.2       httr_1.4.1           backports_1.1.5      assertthat_0.2.1    
##  [26] Matrix_1.2-18        fastmap_1.0.1        lazyeval_0.2.2       cli_2.0.1            later_1.0.0         
##  [31] htmltools_0.4.0      prettyunits_1.1.1    tools_3.6.2          igraph_1.2.4.2       coda_0.19-3         
##  [36] gtable_0.3.0         glue_1.3.1           reshape2_1.4.3       cellranger_1.1.0     vctrs_0.2.3         
##  [41] nlme_3.1-144         crosstalk_1.0.0      xfun_0.12            ps_1.3.2             rvest_0.3.5         
##  [46] mime_0.9             miniUI_0.1.1.1       lifecycle_0.1.0      gtools_3.8.1         zoo_1.8-7           
##  [51] scales_1.1.0         colourpicker_1.0     hms_0.5.3            promises_1.1.0       Brobdingnag_1.2-6   
##  [56] parallel_3.6.2       inline_0.3.15        shinystan_2.5.0      yaml_2.2.1           gridExtra_2.3       
##  [61] loo_2.2.0            StanHeaders_2.21.0-1 stringi_1.4.6        dygraphs_1.1.1.6     filelock_1.0.2      
##  [66] pkgbuild_1.0.6       storr_1.2.1          rlang_0.4.4          pkgconfig_2.0.3      matrixStats_0.55.0  
##  [71] evaluate_0.14        lattice_0.20-38      labeling_0.3         rstantools_2.0.0     htmlwidgets_1.5.1   
##  [76] processx_3.4.2       plyr_1.8.5           magrittr_1.5         R6_2.4.1             generics_0.0.2      
##  [81] base64url_1.4        txtq_0.2.0           DBI_1.1.0            pillar_1.4.3         haven_2.2.0         
##  [86] withr_2.1.2          xts_0.12-0           abind_1.4-5          modelr_0.1.5         crayon_1.3.4        
##  [91] rmarkdown_2.1        progress_1.2.2       readxl_1.3.1         callr_3.4.2          threejs_0.3.3       
##  [96] reprex_0.3.0         digest_0.6.23        xtable_1.8-4         httpuv_1.5.2         stats4_3.6.2        
## [101] munsell_0.5.0        shinyjs_1.1